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We study the nonlinear propagator, a key ingredient in renormalized perturbation theory (RPT) 
' that allows a well-controlled extension of perturbation theory into the nonlinear regime. We show 

that it can be thought as measuring the memory of density and velocity fields to their initial 
5_l ■ conditions. This provides a clean definition of the validity of linear theory, which is shown to 

be much more restricted than usually recognized in the literature. We calculate the nonlinear 
propagator in RPT and compare to measurements in numerical simulations, showing remarkable 
agreement well into the nonlinear regime. We also show that N-body simulations require a rather 
large volume to recover the correct propagator, due to the missing large-scale modes. Our results 
C^l , for the nonlinear propagator provide an essential element to compute the nonlinear power spectrum 

in RPT. 
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How long do perturbations "remember" their initial state? What sets the validity of linear perturbation 
theory? The answer to these questions is important because, for example, measures of the power spectrum in 
the linear regime are being used and planned to be used at higher redshifts @] to determine cosmological 
parameters. In the linear regime the fluctuations being observed contain the full information in the initial 
conditions fl3| . The breakdown of linear theory means that, unless we can somehow model deviations from 
+3 . it, only a very restricted range of scales is reliable for cosmological parameter estimation. 

In this paper we study the propagator of density and velocity fields, the main ingredient that enters into 
a well-controlled extension of perturbation theory into the nonlinear regime (see our companion paper , 
hereafter paper I) . We shall see that the propagator measures the memory of initial conditions and its decay 
at high wavenumbers k is nearly Gaussian with a characteristic scale that provides a clean definition of the 
scale at which linear theory breaks down. The end result of our study is an analytic expression for the fully 
nonlinear propagator of density and velocity fields that can be used in the formalism developed in paper I 
to calculate the nonlinear power spectrum [8j, and other statistics. 

The most common criterion used in the literature to decide the validity of linear perturbation theory is 
to calculate the linear power spectrum and look for the scale at which the rms fluctuations are unity. This 
is intuitively reasonable, but in reality is very unreliable. The reason is that nonlinear corrections depend 
sensitively on the shape of the initial power spectrum (see e.g. Fig. 12 in 4]), not just on its amplitude. 
A better approach is to see how well does the final configuration of density and velocity fields looks like 
evolved from initial conditions by just pure linear theory. One statistical way of doing this is by calculating 
the cross-correlation coefficient r„ between initial and final conditions as a function of scale, 
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r a (k) ee - Xa{k) (1) 
y/Pa(k)P (ky 

where ( 4 r a (k)<5o(k') } = X a (k) <5D(k + k') is the cross-spectrum between initial (described by the density 
So, velocities are assumed in the linear growing mode) and final conditions (denoted by with a = 1 




FIG. 1: Cross-correlation coefficients between initial and final conditions r a as a function of wavenumber k for density 
(a = 1, solid lines) and velocity divergence (a = 2, dashed lines) fields at redshift z = 3 (left panel) and z — (right 
panel). Initial conditions are at redshift z — 35. 



for density and a = 2 for velocity divergence field, see Eq.J2J) for more details). In Eq. Q P$ denotes the 
initial power spectrum, and P a the final one. Figure ^ shows r a {k) at redshift z = 0,3 for density and 
velocity fields. Note that, as expected, the cross-correlation decays slower for z — 3 than z = 0, indicating 
that linear theory is valid over a larger range of scales at z = 3 than z = 0. However, the characteristic 
scale of decay does not differ by as much as one would guess from the condition that the fluctuations be 
of order unity. Indeed, if one estimates the nonlinear scale k n \ from the linear spectrum, iirk^PL(k n i) = 1, 
k n i(z = 0) = 0.21 /iMpc - and k n \(z = 3) = 4.02ft,Mpc - , a ratio of about 20, whereas the characteristic 
scale of decay of cross-correlations between initial and final fields (though it agrees with fc n i at z = 0) only 
shifts by a factor of about 2 — 3 at z — 3. Therefore, the "non-linear scale" can be misestimated by a large 
factor if 47rfc^ 1 Pz ( (fcni) = 1 is used at high redshifts. 

In practice, the cross-correlation coefficient between initial and final conditions is not entirely satisfactory 
as a measure of deviations from linear evolution, since the denominator in Eq. Q is also sensitive to 
nonlinearities through P a (k), and this can lead to cancellations or enhacements. In particular, from Fig.^ 
it appears that the density becomes nonlinear at slightly larger scales (lower k) than the velocities; however, 
this turns out to be incorrect, being more a reflection of the fact that nonlinear corrections to P a (k) are 
mostly positive for the density (a = 1) and negative for the velocities (a = 2). As we shall discuss in detail 
below, this problem can be avoided by the use of the propagator of density and velocity fields rather than 
the cross-correlation coefficient, the propagator G a b is defined by 0,0 

G abM 5 D (k-kO^^^, (2) 
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where ^ a represents the final conditions and <j) a the initial conditions, and again a — 1 corresponds to 
density and a — 2 to velocity divergence fields. As we show below, for Gaussian initial conditions the 
propagator is proportional to the cross-correlation coefficient with a proportionality constant that gets rid 
of the dependence on P a (fc), avoiding the problem with r a . It is obvious from the definition in Eq. (J2J) that 
if linear theory is valid the derivative yields a scale-independent result (the growth factor, for growing-mode 
initial conditions), and any dependence on scale is due to nonlinearities. As k increases nonlinear effects 
become important and the final state does not "remember" the initial one, driving the derivative to zero on 
average. In this sense, the propagator can be though as a measure of the memory that perturbations have 
to their initial conditions. The characteristic scale of decay of the propagator is the cleanest measure of the 
scale where nonlinearities become important. 

This paper is organized as follows. In the next section we review the basics of perturbation theory, see 0,0 
for more details. Section IlIII describes how we arrive at the analytic expression for the nonlinear propagator 
in RPT. Section IIVI presents how to measure the propagator in numerical simulations and compares the 
results with RPT. Section El presents our conclusions. 



II. EQUATIONS OF MOTION AND PERTURBATION THEORY 



In this section we briefly review perturbation theory (PT) in a form useful to facilitate the process of 
resummation (see 7] for more details). The equations of motion for the evolution of dark matter can be 
written in a compact way by introducing the two-component "vector" 

* a (k,ry) = (5(k, V ), -e(k,T,)/H), (3) 

where the index a = 1,2 selects the density or velocity components, with <5(k) being the Fourier transform 
of the density contrast <5(x, r) = p(x)/p — 1 and similarly for the peculiar velocity divergence 6 = V • v. 
Tt = dlna/dr is the conformal expansion rate with a(r) the cosmological scale factor and r the conformal 
time. The time variable ry is defined from the scale factor by 



7] = lno(r), 



(4) 



corresponding to the number of e-folds of expansion. Here we consider a cosmology where VL m = 1 and 
f^A = 0, see section All Fl for the more general case. The equations of motion in Fourier space can then be 
written as (we henceforth use the convention that repeated Fourier arguments are integrated over) 



fi^* (k,» 7 ) + n 6* 6 (k,t ? )=7g J (k ) ki,k 2 ) *„(ki,77) * c (k 2 ,ry), 



where 



-1 
-3/2 1/2 



(5) 



(0) 



(s) 

and the symmetrized vertex matrix 7^ c describes the non linear interactions between different Fourier modes 
and is given by 



7 W(k,ki,k a ) = fo(k-ki-k a ) 

7222( k : k l> k 2) = fo(k-kl-k 2 ) 



(ki + k 2 ) -ki 

|k! +k 2 | 2 (k! -k 2 



(7) 
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7^ c (k, ki, kj) = 7^ b (k, kj,ki), and 7 is zero otherwise, <5d denotes the Dirac delta distribution. The formal 
integral solution to Eq. (JSJ is given by (see [1, IE f° r a detailed derivation) 



rV 

tf „(k, i,) = g ab (r)) 06 (k) + / drf g ab (r] - rf) 7^( k . k i, k 2) *c(ki, »/)*«i(ka, ?/), 

Jo 



(8) 



where r/> a (k) = x J , a (k,?7 = 0) denotes the initial conditions, set when the linear growth factor a(r) = 1 and 
f] = 0. The linear propagator g a b{v) is given by 



9ab(v) = y 



3 2 
3 2 



,-3^7/2 



-2 2 

3 -3 



(9) 



for 77 > 0, whereas g a b(il) = f° r ?7 < due to causality, and g a b{il) — > ^ab as 77 — > + . The linear propagator 
is the Green's function of the linearized version of Eq. |J5J| and describes the standard linear evolution of the 
density and velocity fields from their initial state. 

An explicit expression for \I/ a (k, 77) can be written in the form of a series expansion 



* a (M) = £*W(M), 



n=l 



with 



*<»>(k, Tj) = J &(k - ki... n ) ^ a2 ... a „ (k x , . . . , k„; V ) O1 (ki) . . . Q „ (k, 



(10) 



(11) 



where ki...„ = ki + . . . + k n . Replacing Eq. (|10fl and l|ll|) into Eq. (JSJ determines the recursion relations 

satisfied by the kernels ^Fia\a 2 ...a„- 

In the standard form of Perturbation Theory (see for a review), only the fastest growing mode is 
considered at each order. In contrast with this, Eq. I|ll|) gives the full time dependence, including all 
transients from initial conditions 0, Q ■ This is important since it plays a key role in allowing the process of 
resummation of the non linear propagator. See paper I for a detailed discussion of this point. 

It is worth noticing that so far we have not yet assumed any particular kind of initial conditions. From a 
physical point of view, the most interesting initial configurations are those where #(k, 77 = 0) and 9(k, n = 0) 
are proportional Gaussian random fields, and thus we can write </> a (k) = u a 5o(k), with <5o(k) Gaussian 
distributed. Of particular importance when setting initial conditions in simulations are the growing mode 
initial conditions, for which u = (1,1). Only the case r/> a (k) = ii a r5o(k) will be considered in this paper. 
Replacing this relation in Eq. (|11|) we obtain a simplified version of it 



*t n )(k,»7)= / fo(k-k 1 ... n )^ n )(k 1 ,...,k n ;» 7 )5 (ki)...io(k n ), 



(12) 



where = Taa x ... an The recursion relations for the kernels j^ n > are given in 0. 
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III. CALCULATING THE PROPAGATOR 



A. The Propagator as a Memory of Initial Conditions 

Non linear interactions lead to deviations from the linear evolution described by the first term in the 
r.h.s of Eq. (JSJ. It is possible to interpret this effect as a modification of the linear propagator defined in 
Eq. 10 , this leads to a propagator renormalization. The study of this renormalized quantity, the non linear 
propagator, is the main goal of this paper. 

From its definition, Eq. (J2J, we see that the non linear propagator quantifies the dependence of the density 
and velocity fields on their initial values, in an ensemble average sense, and as we now show it has a very 
simple physical interpretation as a measure of memory of initial conditions. With the help of Eq. (|1(J|) we 
write Eq. as a series expansion 

Gab(k, v) = 9ab{V) + 2_, ( ^ b (k) / ' ^ > 

where we have explicitly separated the linear part from the non-linear contributions. Now we want to make 
use of Eq. 031) for *i n) , assuming Gaussian initial conditions. This implies that the statistical properties of 
the fields </> a (k.) are completely characterized by the two-point correlator 

(0 o (k) b ( k ')> =M k + k ') u a u b P (k), (14) 

where Po(k) denotes the initial power spectrum of density fluctuations. Therefore, for Gaussian initial 
conditions, only the odd terms in the series for G ob will give nonzero contributions after taking the functional 
derivative on vf^™) and performing the ensemble average. Replacing Eq. into Eq. (J2J, we arrive at 

G ab (k,r 1 )=g ab (r,) + J2( 2n + 1 ) l} - / ^ 2 b " +1) (k, k x , -k x , . . . , k„, -k n ; rfiP^) • ■ ■ Po(k n ), (15) 

71 = 1 

with 

•^lr +1) ( k ' k l'- k l' • ■ ■ > k «, -k n ;r/) = ^iba^L,.^' 1 ^' _kl ' ■ ' ■ ' k ™' ~ k «^) U &l ■ ••«ajn' (16) 

As pointed out in the introduction, the non linear propagator is closely related to the cross-correlation 
between final and initial conditions (assumed to be Gaussian). To see this, we multiply Eq. I|llf) by c (k'), 
and ensemble averaging assuming Gaussian initial conditions it follows that 

G ab (k, v ) <<Mk)0 c (k')) = (* a (k, v) <M k ')>- (17) 

Therefore, while in the non linear case it is not true anymore that ty a (k,i]) = G a b(k, r/)06 (k) , since G a b is 
not the Green's function of Eq. JSJ, the nonlinear propagator plays the role of a Green's function in two- 
point sense, (^atyb) — G ac {4> c (j)b}. This is why G ab plays an essential role in calculating different correlation 
functions at the non linear level 0- Also note that Eq. i|17fl makes explicit that G a b cannot depend upon 
the direction of k, since both the initial and final configurations are statistically homogeneous and isotropic. 

For the particular case when </>f,(k) = Uf,5o(\i) and growing mode initial conditions, u = (I, I), it follows 
that 
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FIG. 2: Diagrams for the non linear propagator G(k, rf) up to two loops. 



fo(k + k') [G al (k,r]) + G a2 (k,r])} 
or, in terms of the cross-correlation coefficient, Eq. (JTJ 



Po(k) 



(18) 



G al (k,rj) + G a2 {k,i]) = r a (/c,r) 1 



' Pg(k,T) 

Po(k) 



(19) 



Therefore, the nonlinear propagator quantifies the memory of a density and velocity divergence field to its 
initial conditions as a function of scale and time, a direct measure of the correlation between the final and 
initial configurations. 



B. Diagrammatic Representation 

In order to calculate analytically the terms in the series expansion of Eq. J5J), it is useful to resort to a 
diagrammatic representation (see paper I for more details). For Gaussian initial conditions, only the odd 
terms in the expansion for G a b are non zero, thus, we define 



^)(M)fo(k-k- )s / w ^^ ), (20) 



so that Eq. becomes, 

G ab (k, v ) = g ab ( v ) +^6G%\k,Ti). (21) 
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In principle, an analytic expression for 5G J can be obtained, given the kernel , by integrating out 
Eq. (|15|l . However, this path becomes very cumbersome with increasing n. In order to overcome this difficulty, 
and simplify the process of resummation, we developed in paper I |7J a systematic approach in which the nth 
term in the series expansion of Eq. Ij21|) is represented by a Feynman diagram with n loops. Then, a set of 
rules enabled us to establish a one to one correspondence between diagrams and the corresponding integral 
contributing to 5G^ (k,rj). All the diagrams for the nonlinear propagator up to two loops are shown in 
Fig. |21 Notice that the tree level diagram corresponds to the linear propagator g a b defined in Eq.©. 

A detailed description of the procedure to draw the diagrams, for both ^>a and SG^ , and the rules to 
write down the corresponding integrals is given in paper I, here we only give a brief summary. The symbol (g> 
represents the initial power spectrum u a P( ) (q)u{ ) [see Eq. (|14|) ]. from where two lines emerge carrying opposite 
wavenumbers, say q and — q. Each line finishes in a vertex, and represents a linear propagation g a b{si) from 
the initial conditions (rj — 0) to an arbitrary time of interaction Si at the vertex. At each vertex two modes 
interact, say ki and k2, to give a nonlinear contribution to a third one, k, with k = ki + k2 (conservation 
of momentum). Each vertex in a diagram represents a vertex matrix 7^(k, ki,k2), defined in Eq. J7J. 
The lines between vertices symbolize linear propagations between the corresponding times, e.g. g(si — Sj). 
In addition, each diagram has an outgoing and an incoming line, identifiable with arrows. They represent, 
respectively, a linear propagator g(r/ — si) and g(s2n), where s\ and S2n are the times corresponding to the 
first and last vertices of a general diagram with n loops. Finally, all the intermediate wave vectors running 
through the loops are integrated over as well as all the time variables Sj, each between [0, 77]. 



C. The One-Loop Propagator 

The first non-linear correction in Eq. (|21|l follows simply from applying these rules and is represented by 
diagram 2 in Fig. |3 It reads, 

SG ab( k >V)=^[ ds i [ ds 2 [ d 3 q g ac (i] - Sl ) 7^ } e (k, q, k - q) gdf(si)u f g eg ( Sl - s 2 ) 
Jo Jo J 

7gw(k-q,-q,k) 9hj(s2)uj g tb (s 2 )Po(q) (22) 

where the Su in the definition of 7W in Eq. (JJI) has been integrated out, making explicit the conservation 
of momentum at each vertex. With the help of Eq. J7J) and assuming growing mode initial conditions, is 
possible to write Eq. I|22|) as 

5G^= dsi ds 2 dx 2irq 2 dq — (l ) g( v - Sl ) T g( Sl - s^" 1 g(s 2 )e s ^ P (q) (23) 

Jo Jo J-i Jo q ^ q ' 



with 



r = 



x k/q + x q/k — 2x 2 l — xq/k 
1-xk/q 



(24) 



and x = k • q/fc q. Equation (12;it can be integrated out to give, 



5G^{k 1 a) = Y J frik)a 



(25) 
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where here m takes the values m — {3, 2, 1, 1/2, —1/2, —3/2} and a(r) = exp(?y), Eq. Q. It is possible to 
show that only 4 functions /• ■ are independent, out of the 24 in the set. We take these four to be, 



f= 5 f (3) „ = ^ f (-3/2) /,= 5,(1/2) 

J — gjn ) y— 2-'ii ' 2 11 



i= -f (1) 



(26) 



The normalization is such that in the large-fc limit /, g, h, i 
they all vanish as fc — > 0. Their explicit expressions are 



\k 2 a 2 , with a 2 = | / P (q)d 3 q/q 2 , while 



/(fc) = 

i(fc) = 



1 



504fc 3 a 5 
1 

168fc 3 g 5 
1 

24fc 3 g 5 
-1 



6fc 7 g - 79/s 5 g 3 + 50g 5 fc 3 - 21fcg 7 + -(fc 2 - g 2 ) 3 (2fc 2 + 7q 2 ) In — 4» 

4 |fc + g| 



6fc 7 g - 41fc 5 q 3 + 2fc 3 g 5 - 3fcg 7 + -(fc 2 - g 2 ) 3 (2fc 2 + q 2 ) In 

4 " + 

6fc 7 fl + fc y + 9kq 7 + hk 2 - g 2 ) 2 (2fc 4 + 5fc 2 g 2 + 3g 4 ) In 

4 " Ife + gp 

6fc 7 g + 29fc 5 g 3 - 18fc 3 g 5 + 27kq 7 + ^(fc 2 - g 2 ) 2 (2fc 4 + 9k 2 q 2 + 9g 4 ) In | , 



f (?) A, 
Po(«)d 3 <z, 



|fc-g| 2 



72fc 3 g 5 

Equation i|25|l can now be written in terms of these independent functions in a remarkably symmetric way 



Po(q) d 3 q, 
(27) 



5G$(k,a) = \aa{a) /(fc) - ~/?(a) i(k) - ? 7 (a) fc(fc) + V 3 / 2 <5(a) .g(fc), 
5 5 5 5 

2 2 2 2 

a) = 7«a(a) /(fc) - -/3(a) &(&) + - 7 (a) ft(jfc) - -a" 3/2 6(a) /(fc), 

O 

SG { 2 \\k,a) = La(a) g{k)-^{]{a) i{k) + l 7 (a) i(k)-la-V 2 6{a) g{k), 
5 5 5 5 

6Gg\k,a) = ^aa(a) g(k) - jj/3(a) fc(fc) - | 7 (a) + ^- 3 / 2 5(a) /(fc), 



(28) 



with 



a(a) = a 2 - |a + V 3 / 2 , /3(a) = |a 2 - a + f a" 1 ^ 7 ( a) = H 2 _ a Va + 3-1/2 £( a ) = | a 7/2 _ L + L 
55 55 5 5 55 

The low-fc limit of the one-loop propagator, = g + 5G^ X \ can be computed straightforwardly from 
Eqs. and reading [| 



^ofc (*'°) = .9a&(a) - (ka v ) z a 



2„3 



61/350 61/525 
27/50 9/25 



(29) 



where we kept only the fastest growing mode. Using Eqs. (J5J) and i|29|) we obtain for growing-mode initial 
conditions 



Gs — G\i + G12 ~ a 
Gg = Gii + G22 ~ a 



1 (fcer„a) 2 

210 V ; 

9 

1 ~ W a > 



(30) 
(31) 
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Tl Tl S 2 S, T) S 4 S, S 2 S, 

g(l) g(s 2 -Sj) 



FIG. 3: Tree diagrams for ^(k, 77), up to 4 vertices, that lead to the non linear propagator's resummable subset. 



where we have defined the density propagator Gs and the velocity- divergence propagator Gg for growing mode 
initial conditions. We see that the one-loop corrections define a characteristic scale for which the propagator 
deviates significantly from its linear (fc-independent) value. We shall come back to this important point 
below after we implement the process of resummation. 

D. Resummation in the large-fc limit 

As we discussed in paper I, the non linear propagator plays an essential role in renormalized perturbation 
theory, since it allows a well-defined expansion for the power spectrum and higher-order correlation functions 
into the nonlinear regime. Therefore, one of the main objectives of this paper is to find an adequate analytical 
description for G a b- In order to achieve this, we would need to sum the infinite number of terms defining 
G a b in Eq. (|21() . equivalent to adding up an infinite number of diagrams, shown in Fig. [5] up to two loops. 

Although it is very likely that this resummation cannot be done exactly, we shall see that reasonable 
physical arguments can be put forward to simplify the task enough so that it becomes doable. The first step 
is to realize that since the propagator measures how memory to initial conditions is lost, the next best thing 
to having an exact calculation of it, is to include in the calculation those processes that lead to the slowest 
decay, or memory loss, since they are the ones that establish the dominant behavior. 

Different "processes" first appear at two-loops, where there is more than one diagram contributing. The 
two-loop diagrams in Fig.|21can be organized according to how many interactions (vertices) take place along 
the principal path of the diagram (that connects the beginning to the end without going through initial 
conditions). For example diagrams 3,4,5 have all four interactions along the principal path, diagrams 8,9,10 
have three and diagrams 6,7,11 have only two interactions along the principal path. Diagrams that have 
all interactions along the principal path are special in the sense that the principal path interacts only with 
waves coming directly from the initial conditions in their growing mode. In addition, in the high-A; limit 
each of these interactions results in a wave still in the growing mode, since as we show below in this case 
JabcUbUc tx u a and this keeps repeating itself interaction after interaction along the principal path. For 
the other diagrams, on the contrary, the interactions outside the principal path do not preserve the linear 
growing mode in the high-k limit (since they involve momenta that cannot be compared to k). Therefore, 
when resummed these should be suppressed compared to the resummation of the ones where all interactions 
happen along the principal path, since they should cross-correlate less with initial conditions. One exception 
to this is when the time between final and initial conditions is small enough, in which case there is no time 
for decorrelation and the distinction we made is not significant. 

The dominant subset as defined above can be resummed exactly in the high-fc limit. In general, the subset 
has (2n — 1)!! diagrams with n loops, and it gives rise to diagram 1 (tree level), diagram 2 (one loop) and 
diagrams 3, 4 and 5 (two loops) in Fig. EI The infinite number of diagrams in this subset originate in tree 
diagrams for \E' Q (k, 77), as described in paper I, with a very simple topology. They all start in two Fourier 
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modes that evolve linearly from the initial conditions in growing mode, with u a = (1, 1), interact (represented 
by a vertex), and give rise to a third mode. This mode, in turn, evolves and interacts with a fourth mode 
that had grown linearly until that time, i.e. with no previous interaction. This process continues until the 
final time n, with 2n interactions needed to describe the n-loop diagrams. In Fig. [3] we show the trees with 0, 
2 and 4 vertices that will originate the tree, one-loop and two-loop diagrams respectively. In contrast with the 
rest of the diagrams, the ones belonging to this subset have the simplest ramification possible and thus most 
direct connection to the initial conditions. As discussed above, we expect this subset to be the dominant, 
although we cannot prove it rigorously since we could not sum up the remaining subset and show it is indeed 
subdominant. 

The expression for $0 , relevant to compute SG)^ in Eq. i(2"TJl . can be written recursively as, 

(k', 3) = 2 f ds n g ab (s - S „) 7 ^(k', q n , k + q^...^ (q„, 8 n )¥^ (k + q^.^, s n ) (32) 
Jo 

with k' = k + q n x and 



« (q, s) = g ab (s) 06 (q) = e s u a <5 (q) (33) 

The factor of 2 in Eq. I|32|l is because all the branchings in Fig. [3] are asymmetric if n > 2. This should be 
replaced by 1 for n — 2. Next, we need to take the functional derivative of Eq. I1U2II with respect to </>&(k). 
An important condition, in order to recover the desired diagrams, is that the functional derivative must be 
taken only on any of the two linear modes that precede the very first vertex in Fig. EI I n other words, 
the derivative must be iterated through Eq. (|32|l until n = 2. After using Eq. for *i 2 " +1) , taking the 
derivative as described and using Eq. (|33|l . we arrive at 

j\I/( 2 " +:L) (k' n) f v f S2n f S2 r 

~TT ,\' - = 2 2 " / ds 2n / ds 2n -i... d Sl g(r) - s 2n ) 7 (s) (k + q 2 „... 1 , q 2 „, k + q 2n -i...i) 
o<f>bW Jo Jo Jo L 

g(*2n - s 2 „-i) • ■ ■ g(s 2 - 5i) 7 (s) (k + qi, qi, k) g(ai)] S (q 2n ) . . . 5 (q x ) e '»»+-+« , (34) 

- ab 

- -(s) (s) / 

where 7 is a 2 x 2 matrix defined as jac = l a b c u b, and k = k + q 2ll 1 . In the large-k limit, where the 
external momentum k is much larger than any internal q i; we can expand 7 to leading order, with the help 
of Eq. 0, as 

I \ 1 k 

7ah (k + q„... 1 ,q„,k + q„_ 1 .j) w -— x n 5 ab (35) 

^ Qn 

where x n — k • q n /kq n , and <5 a f, is the 2x2 identity matrix. Equation l|35l) shows that in the high-fc limit 
labcUbUc oc u a at each vertex for this subset of diagrams, as discussed above. A great simplification of 
Eq. (|34|l follows from the fact that 7 is proportional to the identity matrix in the large-k limit. Replacing 
Eq. I|35|) into Eq. (|34|l . and using that g(u — s)g(s — z) = g(u — z), we obtain 

/ 5\E'i 2 ™ +1) (k', rf) \ (a(r))~l) 2n f k 2n 

= gabirj) TTTTt / d3< ?2« • • ■ ^ qi xi n ... x x (<5o(q 2n ) • • ■ <^o(qi)> (36) 



\ 5<pb(k) I (2n)! J q 2n ---qi 

where we have integrated the 2n time variables according to 
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r d S2n n n d S2n - 1 ... r e*«+-+«i = (a( f o j )2 " 

Jo Jo Jo ( 2n V- 

with 77 = ln(a). The ensemble average of the 2n Gaussian fields leads to (2n — 1)!! equal contributions to 
Eq. (|36|l . corresponding to all different pairings of the fields 8q. In addition, for any given pairing there are 
n equal contributions, each given by 

f [ — X iXj (SoidM^)) d 3 q t d 3 qj f ^fd 3 q = -k*o* (37) 

Replacing Eq. ij!T7|) to the nth power into Eq. (|3*S|l . and the result into Eq.J2U we get 

G ab (k, a) = g ab (a) + g ab (a) £ [~k 2 <7 2 v (a - l) 2 ] " (38) 

This series can be summed up to give 

G ab (k, a) = g ab {a) exp(-k 2 a 2 (a _ if/2), (39) 
which shows that we expect the propagator to decay as a Gaussian in the large-fc limit. 

E. The Propagator in Renormalized Perturbation Theory 

We would like now to extend the result in Eq. (|39|l to obtain an analytic result for the nonlinear propagator 
than can be used at any k. The Gaussian dependence on k is also recovered, as shown in paper I, for the 
exact density propagator Gs — G11 + G12 in the Zel'dovich approximation, which has a Gaussian dependence 
on k, for all k, Gf A (fc, a) = ae~? k <T « a for growing-mode initial conditions. 

Based on these considerations we assume that the nonlinear propagator has a nearly Gaussian behavior, 
with a characteristic scale of decay that can be then determined from knowing the behavior of the propagator 
for \ow-k, the one-loop result in Eq. (J2SJ. 

In order to do this and obtain a well defined nonlinear propagator we first rewrite SG^ in Eq. 128(1 in 
a way similar to the linear propagator. That is, we group the contributions to each component of SG^ in 
two terms, one "growing" (with an overall factor a), and one "decaying" (with an overall a~ 3 / 2 ). The only 
guidelines for defining each term are that the most growing contribution (i.e. a) must belong to the growing 
part. In addition, both terms must be always negative, in particular in the limits of large a (for all scales) 
and large-fc (for all times). This guarantees the nonlinear propagator will decay in the long-time or large-fc 
limit. 

From their definition in Eq. J57J), it is easy to check that f,g,i are always negative and diverge 
monotonously to —00 as fc increases, as — fc 2 <7 2 (r)/2. Moreover, h(k) > f(k) > g(k),i(k) for all fc. Provided 
with this and the definitions of a, j3, 7 and 8 in Eq. I|29|l . it is possible to conclude that the only feasible 
way to satisfy our requirements is to group things in the following way, 

511(a) + 5G[\\k, a) = ^a + l a~ 3/2 + | a [a(a) /(fc) - p g (a) i(k) } + \ a^' 2 [5(a) g(k) - ld (a) h(k) ] 

ODD D 
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5n - 


h<5G« - 


-> Gn(k,a 


512 - 


hoG« - 


-> G 12 {k,a 


521 " 




-> G 2 i(fc,a 


522 " 


\-8G$ - 


-> G 22 (fc,a 



512(a) + SG$(fc, a) = I a - I a" 3 / 2 + | a [a(o) /(fc) - /3 9 (a) h(k) ] - | a" 3 / 2 [6(a) f(k) - 7d (a) fe(fe) ] 

ff21 (a) + 5G^ 3 (fc, a) = I a - | a" 3 / 2 + | a [a(o) + 7 9 (o) fc(*0 ] - | a" 372 [6(a) g(k) + f3 d (a) i(k) } 

522(a) + 5Gg\k, a) = I a + | a" 3 / 2 + | a [a(o) <?(&) - \ lg (a) i(k) ] + | a^ 2 [6(a) f(k) - ^/3 d (a) h(k)], 

(40) 

where we have redefined (3 and 7 in Eq. (|29(l through aj3 g (a) = ar z / 2 j3d{a) = /3(a), and aj g (a) = 
a>~ 3 ' 2r fd(a) = 7(a). Next we regard Eq. (f^7Tf> as a nearly Gaussian expansion of the full result and write 

5 ae («(a)/(fe)-/3 9 (aWfe)) + ? a -3/2 e {f(a) S {*)-7 i (a)h(i)) ) 

5 5 

? ae (<»(«)/(fe)-ftW''W) _ 1 ffl -3/2 e («(a)/(k)- 7 «i(a)h(fc)) 

5 5 

^ oe («(«)ff(fc)+7s(o)/i(fc)) _ ^ a -3/2 e (5(a)g(fc)+/3 d (a) l (fc)) ; 

5 5 

H oe (a(a)g(fc)-(3/2) 7B (o)t(k)) + ! fl -3/2 e ({(a)/(k)-(2/3)Wa)M*)), (41) 
o 

This is our prediction for the nonlinear propagator in RPT that we will contrast against N-body simulations 
below. Note that there are no free parameters, given an initial spectrum and cosmological parameters we 
can fully predict G a b(k,rj). Before we do so, however, we must take into account how our calculation is 
changed when we go beyond the £l m = 1, Q\ = model we have considered so far. 

F. Dependence on £l m and Oa 

The dependence on cosmological parameters ignored so far can be taken into account, to a very good 
approximation. First, if we redefine 

¥ (k, rf) = (*(k, ij), -0(k, v )/nf) , (42) 
where / = d In -D+/d In a, with D+ the linear growth factor in the appropriate cosmology, and 

V = In D + (t), (43) 
the equations of motion take the same form as in Eq. (JSJ with the only change 



&ab = 



-1 

-zn m /2p 1/2 



(44) 



These results are so far exact. The problem in Eq. i|44|) is that both Q m and / are functions of time, 
therefore the Laplace transform solution in Eq. |JSj is not correct anymore. Although the growing mode 
will be correct, due to the change in Eqs. (|42I43II . the decaying mode is not, since in general it is given by 

3 / 2 

the Hubble constant H which does not scale as D, ~ in the general case. However it is well known in the 
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standard formulation of perturbation theory, that the dependence of the perturbative solutions in the values 
of fl m and f^A is extremely weak, once the linear growth factors have been scaled out 0, 0, ^| . In other 
words, almost all of the information about the cosmological parameters is encoded in the linear growth factor 
D+(t). The reason is that during most of the time evolution Q m / f 2 ss 1 |9j. Therefore, when we compare 
the propagator in RPT against simulations for the ACDM model in the next section, we will simply replace 
a(r) in Eq. (|41|l by the corresponding linear growth factor D+(r, f2 m , CIa)- 

IV. MEASURING THE PROPAGATOR IN NUMERICAL SIMULATIONS 

We now describe how to measure the propagator in numerical simulations. The simulations we use 
here correspond to three different box sizes, Lbox = 100, 239.5, 479/i~ 1 Mpc with a number of particles 
AT par = 200 3 , 256 3 , 512 3 , respectively. They all correspond to the ACDM model with fl m = 0.3, fl A = 0.7 
and G 8 = 0.9 at redshift z = 0. The small and medium simulations have an input power spectrum with 
shape parameter T — 0.21, the large simulation has instead an input power spectrum that includes the full 
transfer function with £1/, = 0.04 in baryons. We have 24 realizations of the small simulation (see 01 f° r 
a description of these), and only one realization of the medium and large simulations. In the following we 
present results from the large simulation, except for section TlV Dl where we explore the dependence on the 
simulation volume. 

The propagator has never been measured in a numerical simulation before. Here we explore two methods, 
one based on an implementation of the numerical derivative involved in its definition, Eq. the other 
based on the relationship between the propagator and the cross-correlation coefficient, Eq. (I19fl . We show 
that both methods give the same result. 

A. Computing the Functional derivative numerically 

According to Eq. @ , the propagator involves a functional derivative of the final conditions with respect 
to the initial conditions, which formally is given by 

ggoOO = j. vM<Mk) + 6fe(k-kQ]-vM<Mk)] (45) 

In order to measure G a b for a given Fourier mode kj, based on this definition we would need to evolve 
two simulations whose initial configurations of density and velocity fields are equal, except at the value of 
the Fourier coefficient <^,(kj), where they differ by a small amount e (but consistent with Gaussian initial 
conditions). Next, the difference of the final Fourier coefficients x P a (ki) must be divided by the initial 
difference e. This needs to be done several times for each Fourier mode, to obtain an expectation value, and 
repeated for each scale of interest. Computationally, this clearly represents a nearly impossible task. 

Instead, we proceed based on an assumption that resembles the notion of "ergodicity" , namely, that 
all Fourier modes kj, within a bin of wave vectors of magnitude around |k|, can be thought as different 
realizations of the same Fourier mode. This is the same idea one uses in calculating the power spectrum 
from a single realization, here also making use of the fact that G a b depends only on the magnitude of k due 
to statistical homogeneity and isotropy (see section fill A|) . as does the power spectrum. Therefore, given 
one realization of the simulation we compute the difference between \E' a (ki) and the value of the field $„ 
at another Fourier mode kj, within the same bin determined by |k|. Then, we divide the result by their 
difference at the initial conditions, 4>b (kj) — 4>b (kj). Finally, we average over all pairs of modes within each 
bin. 
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k [h Mpc" 1 ] k [h Mpc- 1 ] 



FIG. 4: The density (left panel) and velocity divergence (right panel) propagators at redshift z = from initial 
conditions at z — 5 (corresponding to D+ — 4.68). The symbols represent measurement in numerical simulations, 
from implementation of the functional derivative (crosses) and from the relation to the cross-correlation coefficient 
(circles). The solid lines show the predictions of one- loop PT, Eq. 



All the simulations we use were set up using growing mode initial conditions, where 6 — 9 = So initially. As 
a result of this, it is only possible to measure a "directional functional derivative" {8^ a / S(j)b) = G a i + G a 2, 
which following Eqs. (|3U13 1|> above we shall call the density (G$ = Gn + G12) and velocity divergence 
(G 8 = G 2 i + G 22 ) propagators. 

In summary, given the Fourier coefficients of the final density (a — 1) or velocity divergence fields (a = 2), 
and the initial conditions <5o(k), we estimate the non linear propagators as 



G al (k, V) + G a2 (k, r,) = j^-^ ^ ^ So^TWi) ' ( } 

where the sum runs over pairs of Fourier modes in the bin of magnitude |k|. The result is real because the 
sum involves pairs of Fourier modes with opposite directions. 



B. Computing the Propagator from the cross-correlation 



An alternative approach is to take advantage of the relationship in Eq. I|18|) between the propagator and 
the cross-correlation coefficient. This leads to a straightforward method to measure the density and velocity 
propagators directly from simulations as 
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G«i(M) +G a2 (k, V ) = ^^(k^^JoC-ki), (47) 

Po(fc) AT ^ 

where the sum runs over Fourier modes in the |k| bin. Notice that the two methods proposed to measure 
G a i + G a 2 in Eqs (|46|) and H47[) are intrinsically very different. One is rooted in the definition of the non 
linear propagator as a functional derivative while the other exploits the relation between G a b and the cross- 
correlation coefficient. Nevertheless, Fig. 0] shows that the two methods give essentially the same answer. 
Since the cross-correlation method is much faster (being of order ./V as opposed to A^ 2 ) we shall use that in 
the following. The only limitation of Eq. I|47|l is that it assumes Gaussian initial conditions, however, as long 
as the initial redshift is large enough this should not be a problem, as Fig. 0] shows for initial = 5. 

Figure 01 also includes a comparison with the one- loop propagators from Eq. J5SJ), showing that in the 
low-fc limit the one-loop result does very well in describing the decay of the nonlinear propagators. We 
also see that the one-loop propagator becomes negative (and very quickly large as k increases), which is 
unphysical given its physical interpretation as measuring the memory to initial conditions. This makes clear 
the main point of RPT discussed in paper I: unless the propagator is resummcd to recover the right behavior 
as k — > oo, the loop expansion in standard PT gives rise to large negative contributions which hinder the 
convergence of PT as nonlinear scales are probed. However, from Fig.0]we see that the one-loop corrections 
are not too far off from describing the full decay of the propagator, therefore an approximate resummation 
as described in section Till El stands a good chance of providing a good description. 

C. The Propagator in RPT vs. N-Body Simulations 

Figure [S] shows the comparison between the nonlinear propagator from Eq. (|41|l and measurements in 
the large numerical simulation, for various redshifts z = 0,2,5 (corresponding to linear growth factors 
D + = 28,11.8,5.9) from initial conditions at initial = 35. The left panel shows the density propagator, 
whereas the right panel corresponds to the velocity divergence propagator. From this we see that, despite 
the approximations involved in the resummation process, Eq. I|41l) describes the dependence of the non linear 
propagator with scale and time remarkably well without introducing any free parameters. The reason why 
we can predict the behavior of the nonlinear propagator down to scales in the nonlinear regime can be 
summarized as, 

i) The one-loop correction already describes the onset of the decay of the propagator very well (Fig.0J), 

ii) In the high-A: and long-time limit one expects the subset of diagrams with all interactions along the 
principal path to dominate over the rest. These can be resummed exactly, Eq. (|39|l . leading to Gaussian 
behavior, 

iii) In order to match these two results and thus complete the resummation for all k, we make the assump- 
tion that the one-loop result is a nearly Gaussian expansion of the renormalized propagator. Then, 
as discussed in section fill El there is a unique way to recover the latter from it's one-loop expansion 
satisfying the right asymptotics, leading to Eq. l|4*l"|) . 

D. Dependence of Propagators on the Simulation Volume 

There is one important aspect of the prediction from the RPT propagator that we have ignored so far, 
namely, its dependence on the large-scale modes. This is important because numerical simulations artificially 
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FIG. 5: Predictions for the propagator of density (left panel) and velocity divergence (right panel) from Eq. il 1 1 li 
against measurements in N-body simulations. The three cases correspond (from left to right in each panel) to 
z = 0, 2, 5, with the initial conditions at Zi n i t i a i = 35. The predictions shown in solid lines have no free parameters. 



truncate large-scale modes due to their finite volume. The predictions from Eqs. Ij28(l . Ij39(l . and (|41|l indicate 
that the dependence on the large-scale modes must be significant, since the characteristic scale of decay is 
determined by er" 1 , which is quite sensitive to the smallest wavenumber fcb ox available in a simulation volume. 
In fact, the analytic predictions presented already in all the figures use the exact value fcbox = 27r/Lbox = 
0.013 /iMpc -1 corresponding to the large simulation box. From Eqs. H30I31(I we see that one can define the 
characteristic scale of decay from the low-fc expansion of the density and velocity propagators as 



k s (z) = v/105/61 ex" 1 {z), (48) 
kg(z) = y/E/9 a^z), (49) 

a 2 v (z) = {4ir/3)D 2 + (z) / P (k)dk, (50) 



which is roughly the scale at which G$/D + and Gg/D + have respectively dropped to e -1 / 2 « 0.6. 

Figure shows the dependence of the propagators on the simulation volume, confirming that there is 
a significant dependence, in agreement with the estimates from Eqs. I|48I49|I . A larger volume leads to a 
stronger decay due to the fact that more modes are present and the increased interactions lead to faster loss 
of memory of the initial conditions. A similar situation holds for the cross-correlation coefficient between 
initial and final conditions presented in Fig. ^ I n other words, the larger the volume allows interactions 
among more modes which leads to a larger scatter of the final densities and velocities about the linearly 
extrapolated initial conditions, which suppresses r a . 
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FIG. 6: The dependence of the measured density (left) and velocity divergence (right) propagators on the volume 
of the simulation used. The symbols denote measurements of the propagator at z = from initial conditions at 
^initial = 5 for different simulation volumes, Lb ox = 100 ft -1 Mpc (pentagons), Lb ox = 239.5 h~ x Mpc (squares), 
L box = 479 h~ x Mpc (triangles). 



Why is such a strong dependence on the simulation volume not seen in the power spectrum measurements? 
The reason is that in a finite volume the propagator decays slower, leading to a fictitiously wider validity of 
linear perturbation theory, but the truncation of low-fc modes also suppresses the mode-coupling terms that 
should enhance the power, therefore these two effects compensate each other. However, that means that the 
dependence of the power spectrum on cosmological parameters will be incorrect at the transition scale, since 
these two contributions have a different dependence on them (their ratio being proportional to the linear 
power spectrum) and therefore cannot completely cancel. 

The importance of this finite-volume effect is stronger for steep spectra, in fact one should correct the 
Lb ox = 479/iMpc~ 1 measurements in Fig. by a shift of about 4% to the right according to Eqs. Q48I - 
150(1 to take into account the difference between the initial spectrum of this simulation as opposed to the 
Lhox = 100, 239.5 ft- -1 Mpc ones, but this does not alter the conclusions significantly. 

Moreover, using Eqs. (|48I50I) one concludes that the decay length of the propagators in the Lbox = 
479/iMpc -1 case is still off by about 4% compared to the infinite volume case that will show a stronger 
decay. An important lesson from this is that simulations designed to study the high-redshift Universe that 
generally use smaller boxes can misestimate nonlinear effects due to lack of large-scale modes very easily, as 
the requirement to get the decay scale of the propagator correctly is independent of redshift and, for example, 
to get kg, kg accurate to 10% requires a simulation of about Lbox = 240 hr 1 Mpc, which can be challenging 
for high-redshift numerical studies. 
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FIG. 7: Comparison between the nonlinear scale defined from the linear power spectrum k n i(z) (short-dashed lines) 
and those derived from the characteristic scale of decay of the nonlinear propagator ks(z) (solid lines), ko (z) (long- 
dashed lines) given by Eqs. (1481491 . At high redshift k n i(z) significantly underestimates the importance of nonlinear- 
ities. 



V. CONCLUSIONS 

In paper I jj we developed renormalized perturbation theory (RPT), in which nonlinear evolution in 
the growth of large-scale structure is represented by Feynman diagrams constructed in terms of three ob- 
jects: the initial conditions (e.g. perturbation spectrum), the vertex (describing non-linearities) and the 
propagator (describing linear evolution). We showed that loop corrections to the linear power spectrum 
organize themselves into two classes of diagrams: one corresponding to mode-coupling effects, the other to 
a renormalization of the propagator. The nonlinear propagator that results quantifies the deviation from 
linear evolution of individual Fourier modes, and allows a well-defined perturbation theory that can probe 
the nonlinear regime. 

In this paper we studied the propagator in detail and showed that its decay into the nonlinear regime can 
be thought as measuring the "memory" of perturbations to their initial conditions, being directly propor- 
tional to the cross-correlation coefficient between final and initial configurations. We developed an analytical 
description of the nonlinear propagator by summing its perturbation series using physically motivated ap- 
proximations. The results are in remarkable agreement with measurements in numerical simulations all the 
way into the nonlinear regime. We also described how to measure the propagator in numerical simulations 
using two very different algorithms that gave consistent answers. The analytic results show that the propaga- 
tor is rather sensitive to the finite volume of simulations, and this is indeed in agreement with measurements 
in simulations of different volume. The requirement of getting the right dependence of the propagator on 
scale turns out to be a stringent one for N-body simulations, particularly at high redshift where typically 
smaller boxes are used. 

Coming back to the questions raised in the introduction about the validity of linear perturbation theory, 
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we can give a quantitative answer for a robust definition of the nonlinear scale by looking at the characteristic 
scale of decay of the propagator of density and velocity divergence fields ks, kg given by Eqs. I|48I49[) . Figured 
shows the evolution of ks, kg as a function of redshift and compared to the nonlinear scale k n \ defined from 
the linear power spectrum, 47rA^ 1 i- > £(fcni) = 1. Although the different definitions agree for the density field 
at z = 0, at high redshift the nonlinear scale k n \ typically used in the literature can be more than an order 
of magnitude off from kg, kg, which we argued here is the cleanest definition of at what scale nonlinearities 
become important. 

Since the non linear propagator plays an essential role in the calculation of correlations functions, the next 
obvious step is to use the results presented here to calculate the non linear power spectrum. In the framework 
of RPT, once the nonlinear propagator is known, the power spectrum is given as a mode-coupling series 0- 
Although in principle many terms are needed in order to cover a large range of scales, the study of how 
baryon wiggles are affected by nonlinear evolution requires only a few terms || • Already, the dependence of 
the propagator on scale found here provides a simple estimate of the scale at which baryon wiggles should be 
significantly suppressed (remaining only as e^ 1 / 2 « 0.6 of their linear amplitude over a smooth component), 
given by k$(z),kg(z) in Fig. [7| 
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